% Process the eej strength derived from satellite data
% and compare it with ground deribed eej indices
%start date 30 OCT 2005
%latest date 10 SEP 2006
% Just copy of eej_champ6, but with not figure outputs - autmated writing
% to a file
clear;
more off;
load c:\manoj\projects\eej\OBS_array.mat OBS
%OBSdH1 description 
%The file OBS is Structured array with stations as listed below
%[1-ABG,2-ETT,3-HYB,4-TIR,5-HUA,6-FUQ,7-MBO,8-GUI,9-AAE,10-QSB,11-GUA,
%12-CBI,13-PND,14-ELT,15-AAE(new),16-NGP,17-UJJ,18-KOD,19_YAP,20-BIK]

% Note please update c:\manoj\projects\eej\used_obs.xls if more data are
% added to OBS_array.mat

%-------- correlate the CHAMP derived EEJ strength and ground delta H
% with LT correction
%--------------------------------------------------------------------
fidd = fopen('c:\manoj\temp\eej_slope.txt','at');

[a,b] = xlsread('c:\manoj\projects\eej\used_obs.xls');

bin5 = [-60 -55 -50 -45 -40 -35 -30 -25 -20 -15 -10 -5 0 5 10 15 20 25 30 35 40 45 50 55 60];
bin10 = [-60 -50 -40 -30 -20 -10 0 10 20 30 40 50 60]; 
bin20 = [-60 -40 -20 0 20 40 60];
% obs_pairs = [2,1;2,3;2,13;2,16;2,17;2,18;4,1;4,2;4,3;4,13;4,16;4,17;4,18;5,6;7,8;9,10;9,14;15,10;15,14;...
%     18,1;18,3;18,13;18,16;18,17;13,1;13,3;19,20];
obs_pairs = [2,3;4,1;5,6;15,14;19,20;18,3];
%obs_pairs = [2,3;2,1];
%Control Center
%[1-ABG,2-ETT,3-HYB,4-TIR,5-HUA,6-FUQ,7-MBO,8-GUI,9-AAE,10-QSB,11-GUA,
%12-CBI,13-PND,14-ELT,15-AAE(new),16-NGP,17-UJJ,18-KOD,19-YAP, 20-BIK ]
eej_data_set =1;%1-old, 2-new
satdata  = 2; % %Pear Current Density = 1; av cur distan 2;
obslt    = [9.9 13.1]; %NOTE the LT are represented as for eg 8.99999 or 9.001 the filter
incrlat  = 10;
bin_x   = bin10;
mindate  = tmjd(2000,01,01,0,0,0,1); % starting date
maxdate  = tmjd(2002,12,31,0,0,0,1)+1; % ending date
sqcor    = 0; %SQ correction
ltcor    = 0; %LT correction
robft    = 1; %Robust outlier elimination
plotopt  = 'r-'; %[b g r k]
plotcol  = 'r';
plotchar = 'r';
ncoeff   = 9;%same for hourly means 20 for minute means
kp_thrashold = 2;
dummy_p = 1;

for ikk = 1:length(obs_pairs),

eq_obs = obs_pairs(ikk,1);
sq_obs = obs_pairs(ikk,2);
data1  = [];    

switch eej_data_set,
    case 1,
        data = load('c:\manoj\projects\eej\eej_par_all.txt');
            switch satdata,
                case 1,
                eej_index = 9;
                case 2,
                eej_index = 10;
            end;
    case 2,
        load c:\manoj\projects\eej\EEJ_Stefan data;
        sqcor    = 0;
        switch satdata,
            case 1,
            eej_index = 5;
            case 2,
            eej_index = 6;
       end;
end;
%INFORMATION
%The EEJ parameters in the file
%1-mjd,2-longitude of eq-crossinfg, 3-zero of current position-south
%4 -north,5-half-width of current - south,6-north
%7minimum current south,8-north,9-pear current density 10-average
%difference 11 - diff to south minumum 12 - diff to north minum, 13-total
%east current ,14 total return current., 15 -f,16 - eq angle, 17-
%eq-direction, 18 - mgnteic equator latitude

    
%load c:\manoj\projects\eej\guacm4.mat index sq_est
        
        deltah = OBS(eq_obs).data - OBS(sq_obs).data;
        k = nanmean(reshape(deltah,[24,length(deltah)/24]),2)';  
%         k=circshift(k',[0,round(OBS(eq_obs).long/15)]); % 
        obs_index = 2;
        data1(:,1) = OBS(eq_obs).mjd;
        data1(:,2) = OBS(eq_obs).data - OBS(sq_obs).data;
%       load c:\manoj\projects\eej\ettcm4.mat index sq_est;
        xtime = 0.5:1:23.5;
% 		xtime = 1:24;
        obslon = OBS(eq_obs).long;
        if obslon > 180,
            obslon = -(360-obslon);
        end;
        st = obslon - (incrlat*6+incrlat/2);
        en = obslon + (incrlat*6+incrlat/2);

[p,s] = polyfit(xtime,k,ncoeff);%In case of minute data change the degree to 20
%----------------------------

ndd = zeros([1,25]);
nloop = 1;
xdata = [];
x1data = bin_x;
nbindata = [];
ceff = [];
cerr = [];
nbin = 1;


for ki = st:incrlat:en,
L = data(:,2) > ki & data(:,2) < ki+incrlat & data(:,1) < maxdate  &  data(:,1) > mindate; %The last filter is to select only data before Dec 31, 2002
IndiaE = data(L,:);
ndata = 1;
ind1 = [];
ind2 = [];
nnodata=1;
for i = 1:length(IndiaE(:,1)),
    ind = findnearest(IndiaE(i,1),data1(:,1));
    [lt1,dummy1] = champ_lt(IndiaE(i,1),IndiaE(i,2),obslon);
    [dummy2,lt2] = champ_lt(data1(ind,1),IndiaE(i,2),obslon); % This is because observatory FDAY ~= CHAMP FDAY
    
    if eej_data_set == 1,
        
         if lt2 >= obslt(1) & lt2 <= obslt(2) & abs(data1(ind,obs_index)) < 300 & data1(ind,1) < max(data(:,1)),% & kp(ind) <= 2.5, %The last fliter take care the nodata (9999) points
            [ut_obs] = (data1(ind,1)-floor(data1(ind,1)))*24;
            [ut_sat] = (IndiaE(i,1)-floor(IndiaE(i,1)))*24;
             ind1(ndata) = ind;
            dummy_mjd(dummy_p) = IndiaE(i,1);
            dummy_p = dummy_p+1;
            f_ch(ndata) = polyval(p,ut_sat); %Getting the expexted EEJ strength from poly coefficients for champ 
            f_ob(ndata) = polyval(p,ut_obs); %Getting the expexted EEJ strength from poly coefficients for obs 
%             f_ch(ndata) = interp1(1:24,k,lt1,'spline'); %Found spline interpolation working better 18.1.06 !
%             f_ob(ndata) = interp1(1:24,k,lt2,'spline');
%             fprintf('Lon = %5.1f i = %d f_ch = %f f_ob = %f lt1 = %f lt2 = %f kp = %3.1f\n', ki,i,f_ch(ndata),f_ob(ndata),lt1,lt2,data1(ind,end));
            ndata = ndata+1;
        else,
            ind2(nnodata) = i;
            nnodata = nnodata + 1;
        end;
    elseif eej_data_set == 2,
            kp1 = IndiaE(i,10);
            cmax = IndiaE(i,5);
            if lt2 >= obslt(1) & lt2 <= obslt(2) & lt1 >= obslt(1) & lt1 <= obslt(2) ...
            & abs(data1(ind,obs_index)) < 300 & data1(ind,1) < max(data(:,1)) ...
            & kp1 <= kp_thrashold & cmax >= 0.03, %The second last fliter take care the nodata (9999) points
            [ut_obs] = (data1(ind,1)-floor(data1(ind,1)))*24;
            [ut_sat] = (IndiaE(i,1)-floor(IndiaE(i,1)))*24;
            dummy_mjd(dummy_p) = IndiaE(i,1);
            dummy_p = dummy_p+1;
            ind1(ndata) = ind;
            f_ch(ndata) = polyval(p,ut_sat); %Getting the expexted EEJ strength from poly coefficients for champ LT
            f_ob(ndata) = polyval(p,ut_obs); %Getting the expexted EEJ strength from poly coefficients for obs LT
%             f_ch(ndata) = interp1(1:24,k,lt1,'spline'); %Found spline interpolation working better 18.1.06 !
%             f_ob(ndata) = interp1(1:24,k,lt2,'spline');
%             fprintf('Lon = %5.1f i = %d f_ch = %f f_ob = %f lt1 = %f lt2 = %f kp = %3.1f\n', ki,i,f_ch(ndata),f_ob(ndata),lt1,lt2,data1(ind,end));
             ndata = ndata+1;
             else,
            ind2(nnodata) = i;
            nnodata = nnodata + 1;
            end;
    end;
end;

if length(ind2) > 0,
    IndiaE(ind2,:) = [];
end;

if length(ind1) > 1,
    obsE = data1(ind1,:);
    X = IndiaE(:,eej_index);
    if ltcor == 1,
        Y = obsE(:,obs_index).*(f_ch./f_ob)';
    elseif ltcor == 0,
        Y = obsE(:,obs_index)';
    end;
    
    if robft == 1,    
%    Robust outlier rejection  -- START
     if length(X) > 3,
     [rs,rstat] = robustfit(X,Y);
     residuals = abs(rstat.resid);
     LLL = residuals == max(residuals); %Just remove the most noisy point
     X(LLL) = [];
     Y(LLL) = [];
     end;
    % Robust outlier rejection  --- END
    end;

    [cof,err] = corrcoef(X,Y);  % correlation & Correction for observatory data 
    ceff(nloop) = cof(2,1);
    cerr(nloop) = err(2,1);
    xdata(nloop) = ki+incrlat/2;
    ndd(nloop) = ndata-1;
    nloop = nloop +1;

    %----------- PLOT SCATTER
if x1data(nbin) == 0,
[c,sfit] = polyfit(IndiaE(:,eej_index),obsE(:,obs_index).*(f_ch./f_ob)',1);
end;


%-----------------END
else,
    x1data(nbin) = NaN;
     
end;
nbin = nbin + 1;
clear ind f_ch f_ob ind1 ind2 ndata nnodata IndiaE obsE ltt1 ltt2 ltt3 ltt4 obst chpt X Y;
end;

x1data(isnan(x1data)) = [];
ndd(length(x1data)+1:end) = [];
ss= num2str(ndd');
fprintf('\n');

fprintf(fidd,'%s-%s %6.2f %6.2f ', OBS(obs_pairs(ikk,1)).code,OBS(obs_pairs(ikk,2)).code,c(1),c(2));
fprintf('%s-%s %6.2f %6.2f ', OBS(obs_pairs(ikk,1)).code,OBS(obs_pairs(ikk,2)).code,c(1),c(2));
fprintf(fidd,'%3d ',length(x1data));
L = x1data == 0;
fprintf(fidd,'%3d %3d %4.2f %10.4e ',x1data(L),ndd(L),ceff(L),cerr(L));
fprintf('%3d %3d %4.2f %e ',x1data(L),ndd(L),ceff(L),cerr(L));

% fprintf('%3d ',length(x1data));
% for ipk = 1:length(x1data),
%     fprintf(fidd,'%3d %3d %4.2f %e ',x1data(ipk),ndd(ipk),ceff(ipk),cerr(ipk));
%     fprintf('%3d %3d %4.2f %e ',x1data(ipk),ndd(ipk),ceff(ipk),cerr(ipk));
% end;
fprintf(fidd,'\n');
fprintf('\n');

end;
more on;